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Abstract. 

We evaluate the average expansion rate of a universe which contains a realistic evolving 
ensemble of non- linear structures. We use the peak model of structure formation to 
obtain the number density of structures, and take the individual structures to be 
spherical. The expansion rate increases relative to the FRW value on a timescale of 
10-100 billion years, because the universe becomes dominated by fast-expanding voids. 
However, the increase is not rapid enough to correspond to acceleration. We discuss 
how to improve our treatment. We also consider various qualitative issues related to 
backreaction. 



PACS numbers: 04.4G.Nr, 95.36. +x, 98.80.-k, 98.8G.Jk 



Evaluating backreaction with the peak model of structure formation 1 
Contents 

1 Introduction 1 

2 Clumpy spacetimes 4 

2.1 Light propagation in a clumpy spacetime 4 

2.2 The expansion rate in a clumpy spacetime 7 

3 The Buchert formalism 9 

3.1 The local equations 9 

3.2 The average equations 11 

3.3 First integrals of the Buchert equations 13 

4 Backreaction with the peak model 15 

4.1 Setting up the model 15 

4.2 Preferred time in structure formation 21 

4.3 The expansion rate 23 

4.4 Improving the model 26 

5 Discussion 28 

5.1 Observational issues 28 

5.2 Non-Newtonian effects 32 

5.3 Conclusion 35 



1. Introduction 

Observations and acceleration. There is overwhelming observational evidence against 
homogeneous and isotropic models of the universe based on standard general relativity 
(i.e. the Einstein-Hilbert action in four dimensions) with ordinary matter (i.e. baryons, 
dark matter, neutrinos and photons). Observations are usually interpreted keeping to 
the linearly perturbed homogeneous and isotropic Friedmann-Robertson- Walker (FRW) 
models, in which case it is necessary to modify gravity or introduce matter with negative 
pressure. The ACDM model, which involves the simplest modification of gravity, the 
cosmological constant, or the equivalent form of matter, vacuum energy, is in agreement 
with various observations of cosmological distances. In FRW models, distance has a 
simple correspondence with the expansion rate (and spatial curvature). 

However, linearly perturbed FRW models do not describe the non-linear structures 
present in the real universe. It is also important to recognise that the only evidence for 
modified gravity or exotic matter comes from observations of cosmological distances, 
interpreted in terms of the expansion rate. For example, deviations from general 
relativity have not been observed in the solar system [1]. (The Pioneer anomaly is 
a possible exception. However, it does not agree with the predictions of the modified 
gravity models which have been developed to explain the cosmological observations.) 
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The situation is quite different from that of dark matter, for which there are 
several independent hues of evidence, such as the motions of stars in galaxies, the 
motions of clusters, the peak structure of the cosmic microwave background (CMB), 
the early formation of structures, gravitational lensing, as well as direct measures of the 
matter density combined with the baryon density given by Big Bang Nucleosynthesis. 
For this reason, constructing alternatives to dark matter requires resort to baroque 
models [2], if it is possible at all. In contrast, in order to explain the observations 
without modified gravity or exotic matter, it is only necessary to change the distance 
scale (or the expansion rate) as a function of redshift. 

The measured quantities from which the expansion rate is inferred can be divided 
into background quantities and perturbations, though measurements of background 
quantities also involve perturbations, apart from the luminosities of Type la supernovae 
(SNe la). (In realistic models with non-linear structures, there is no simple division 
into background and perturbations.) Most of the information on the expansion rate 
comes from measures of the background geometry. While the expansion rate has been 
measured from different physical systems, such as the CMB [3], large scale structure [4] 
and SNe la [5], they all probe essentially the same distance measure at different redshifts. 
(The luminosity distance and the angular diameter distance are directly related to each 
other. In FRW models, they can be expressed in terms of the proper distance and 
redshift. For the different distance measures in the context of FRW models, see [6,7].) 

The details of the expansion rate (or the distance scale) are not well known except in 
specific models. Model independent constraints are relatively weak. The CMB is mostly 
sensitive to the angular diameter distance to the last scattering surface, which is related 
to the position of the peaks in the angular power spectrum [8,9]. Type la supernovae 
provide a more direct measure of the expansion rate, but present data is not sufficiently 
accurate to give detailed information. Quoted constraints on, for example, the equation 
of state are often driven by the assumed parametrisation, which can be responsible for 
misleading apparent precision and artificially small confidence level contours. For the 
importance of parametrisation, see [7,10-14]; an analysis of the current situation in 
terms of a piecewise constant equation of state is given in [15]. Additional systematic 
effects, such as metallicity [16], changes in the treatment of dust [17] or the light- 
curve fitting method [18] can further degrade the reliability of the SN la measurements. 
It is noteworthy that the quality of fit of the ACDM model has decreased with the 
introduction of each new SN la dataset up to and including the 'gold' sample and the 
ESSENCE data [19]. This may hint at inadequacy of the ACDM description, or it can 
be indicative of underestimated systematic errors. 

Observational constraints from perturbations are much weaker than those from 
background quantities. The Integrated Sachs- Wolfe (ISW) effect has been detected 
via cross-correlation of the CMB and matter tracers at different redshifts [20]. This 
is interpreted as time evolution of the gravitational potential, which may be due to 
accelerating expansion or spatial curvature. The amplitude is slightly larger than 
expected in the ACDM model, though not significantly so (around 2a). Other 
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constraints from the evolution of perturbations are rather weak [21]. A measure which 
combines background evolution and the evolution of perturbations is the number count 
of clusters as a function of redshift: in an accelerating model the number density should 
rise sharply with increasing redshift. It has been earlier argued that this is not seen in 
the data, and that the observations instead prefer deceleration [22]. However, it seems 
that limited understanding of gas physics prevents drawing reliable conclusions from 
cluster counts at present [23]. 

In the context of homogeneous and isotropic models, we can with confidence 
say that the expansion has accelerated within the last few billion years (though it is 
not necessarily accelerating at present), and determine the overall magnitude of the 
acceleration [11-13, 18]. In order for a model to agree with the observations, it is not 
necessary to reproduce the expansion history of the ACDM model in detail, only to 
produce a roughly similar amount of acceleration in the same era. 

Structure formation and the coincidence problem. While there is no evidence apart 
from the expansion rate for modified gravity or negative pressure matter, we know that 
the assumption of linearly perturbed homogeneity and isotropy breaks down in the 
universe at late times. One has to evaluate the effect of non-linear structure formation 
on the expansion rate (or the distance scale) before concluding that it is necessary to 
change either gravity or the matter content. 

It was suggested in [24,25] that inhomogeneities related to structure formation could 
be responsible for accelerated expansion (the possibility had been earlier touched upon 
in [26,27]). This was discussed more concretely in [28,29] where it was demonstrated with 
a toy model how the formation of non-linear structures can lead to average acceleration, 
even when the expansion locally decelerates (see also [30]). The physical reason is simply 
that the fraction of volume in faster expanding regions rises, so the average expansion 
rate can rise. The bigger is the difference between the slower and faster expanding 
regions, the more rapid is the change as the faster regions take over. It was pointed out 
in [28, 29] that the growth of non-linear structures involves a growing variance in the 
expansion rate, and that structure formation has a preferred time around 10-100 billion 
years, near the observed acceleration era. The timescale emerges from the CDM power 
spectrum, essentially from the time of matter-radiation equality encoded in the change 
of slope of the CDM transfer function. This could solve the coincidence problem of 
why the acceleration has started recently in cosmic history, something that the ACDM 
model does not explain. 

The link between the preferred time in structure formation and the change of the 
expansion rate was only conjectured in [28,29]. The toy model with acceleration involved 
only two regions instead of a realistic ensemble of structures, and had no link to the 
preferred time. In the present work, we remedy both shortcomings. We first discuss 
the effect of structures on the propagation of light and the expansion rate in section 2, 
and set up the Buchert backreaction formalism in section 3. In section 4 we present 
our model for a statistically homogeneous and isotropic universe containing an evolving 
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ensemble of non-linear structures. We calculate the average expansion rate and find that 
it grows relative to the FRW value around 10-100 billion years, though the change is not 
rapid enough to correspond to acceleration. We explain the reason and consider how 
improve our treatment. In section 5 we discuss some observational and theoretical issues 
related to backreaction, and summarise our results. This work is a follow-up to [29], and 
more discussion of these topics can be found there. However, we have repeated some 
material to make the presentation self-contained. 

2. Clumpy spacetimes 

2.1. Light propagation in a clumpy spacetime 

The overall geometry. Though the evidence for negative pressure matter or modified 
gravity is often phrased in terms of the expansion rate, observations of light just 
provide constraints on the distance scale at different redshifts. The conclusion that 
the observations imply accelerated expansion has been established only in the context 
of linearly perturbed FRW models, where the distance scale is directly related to the 
expansion rate (and spatial curvature). The real universe is not locally perturbatively 
near homogeneity and isotropy at late times, but contains non-linear structures (we 
refer to such a universe as clumpy), so the FRW results for light propagation cannot be 
straightforwardly applied. 

Since non-linear structures affect the propagation of light, it might be possible 
to explain the observations without recourse to accelerated expansion. We will be 
interested in the possibility that non-linear structures lead to actual acceleration. But 
even in this case, one needs to study light propagation in a clumpy space to be able to 
compare the model to observations. 

In FRW models, in order to convert between measures of distance and the expansion 
rate, one needs to know the spatial curvature. From the volume and spatial curvature 
of each hypersurface of proper time, the distances can be reconstructed. In a general 
spacetime, this is not true, and the distance scale does not necessarily correspond to the 
expansion rate (see e.g. [31]). We can divide the effect of structures on the passage of 
light into an overall part, which can be expressed in terms of the average geometry given 
by the scale factor (which measures the volume of the spatial hypersurface) and the trace 
of the spatial Ricci tensor (which is a measure of the curvature of the hypersurface), 
and a part which cannot be expressed in terms of the average geometry. 

For example, in perturbed FRW models, the ISW and Rees-Sciama effects are not 
expressible in terms of the scale factor and spatial curvature. These effects are small 
both because the perturbations are small and because there are cancellations present. In 
particular, in the linearly perturbed Einstein-de Sitter universe (the spatially flat FRW 
dust model), the ISW effect is zero independent of the magnitude of the gravitational 
potential (as long as the second order contribution can be neglected). In models which 
include realistic, non-perturbative structures, there is no obvious reason for such effects 
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to be small. Nevertheless, it is possible that this is the case in a universe which is 
statistically homogeneous and isotropic, and that the passage of light can (up to small 
corrections) be described in terms of the overall geometry. 

Studies of light propagation. There is currently no derivation of light propagation for 
the case of a clumpy universe that contains realistic evolving structures like those that 
are actually present. However, several models with realistic elements have been studied. 
The effect of inhomogeneities on the passage of light was first discussed by Zel'dovich 
in 1964 [32] (early literature on the topic also refers to an unpublished colloquium by 
Feynman in the same year). This and the work which followed [33-35] studied the 
passage of light in the case when the deviation from the FRW models is in some sense 
small (see [36-42] for later studies of perturbed FRW universes). 

The effect of non-linear perturbations on the redshift and the luminosity distance 
was studied in 1969 in the context of the Swiss cheese model [43], where sections of 
FRW universe are cut out and replaced with the Schwarzschild solution of equivalent 
density. It was found that corrections to the FRW results could be sizeable. 

The approach introduced by Zel'dovich in [32], where light rays encounter only 
a fixed fraction of the mass of the universe because of clumping, was used by Dyer 
and Roeder in 1972 to calculate corrections to the luminosity distance [44] (and is 
now known as the Dyer- Roeder formalism). They too concluded that large effects are 
possible. The formalism was generalised to include a mass fraction of clumps which 
varies with redshift in [45]. For applications of the Dyer- Roeder formalism and Swiss 
cheese models to observations of SNe la, see [46,47]. 

It was argued by Weinberg in 1976 that the effects of dumpiness on the luminosity 
distance cancel, and the FRW formula can be used to describe light propagation [48]. 
The main argument is that the number of photons is conserved and gravitational 
deflection conserves photon energy, so the luminosity distance can be determined simply 
in terms of the area of a sphere drawn around the emission point. However, it was 
pointed out in [49] that this argument assumes that the area is the same as in the FRW 
model, which is, in fact, the issue under study. An exact counter-example was provided 
in [50], and further arguments that the FRW result may not apply were given in [51]. 

In addition to purely analytical work [52-55], light propagation has been studied 
using ray-tracing methods [38,56-60]. Ray-tracing studies disagree on the magnitude of 
the effect on the passage of light, presumably due to different modelling assumptions. 
Among the most realistic studies are modified Swiss cheese models where the non-linear 
density concentrations reside in the walls around the holes, instead of the center. Studies 
of realistically sized structures have found the corrections to be small, at the percent 
level [58,59]. In [58] the integrated effect of an ensemble of voids was found to be 
proportional to void size divided by the horizon size. A sizeable correction was found 
in [60], but this is consistent with the previous results, since the structures studied in [60] 
are unrealistically large. However, the amplitude of the effect has not been conclusively 
settled, since none of the models have been completely realistic, in particular with regard 
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to the time evolution of structures and spatial curvature. 

In analytical work, there is also no agreement as to the whether the effect can 
be large. It was found in [54] that for an ensemble of uncorrelated, static clumps of 
matter in a FRW universe, the corrections are small (see [53, 55] for discussion of the 
small effects which could be observable). According to [49,51], large effects are possible 
due to the formation of caustics; see also [52]. (For studies along related lines, see the 
program of observational cosmology [61].) 

As summarised in [62] (see also [63]), one important issue is that light propagation 
is in general affected both by the Ricci tensor given by the local matter distribution and 
the non-locally determined Weyl tensor. In FRW models, the Weyl tensor is zero and the 
geometry of spacetime is entirely determined in terms of the local matter distribution. 
In contrast, in the real universe, light mostly travels in vacuum where the Ricci tensor 
is zero, and the geometry is determined by the Weyl tensor. This difference is related 
to the role of shear of the null geodesies, generated by nearby matter. In this sense, the 
propagation of light in the real universe is the reverse of the FRW situation. (See [34,64] 
for early discussion.) 

In contrast to the studies of an ensemble of structures which we know to exist in 
the universe, there has recently been revived interest in the idea that we would live near 
the centre of a single untypically large spherical void ('Hubble bubble'). Studies of the 
passage of light in the spherically symmetric Lemaitre-Tolman-Bondi (LTB) model [65] 
(see [66] for a review) show that if the local bubble is sufficiently large, its effect on the 
passage of light could explain the observations of SNe la (and possibly be consistent 
with the location of the CMB peaks). For references on the 'Hubble bubble' and the 
passage of light in the LTB model, see [29] ; for a review on explaining the observations 
with the LTB model, see [67]. 

The results for the passage of light are not conclusive, and no models have properly 
dealt with the kind of nested, evolving structures that are a central feature of hierarchical 
structure formation. The effects which have been found for realistic models of randomly 
distributed small structures are small. Large effects have been demonstrated only when 
the observer occupies a special location, such as near the center of a large spherical 
region. It therefore seems plausible that those effects of dumpiness on the passage 
of light which have been studied are small if the matter distribution is statistically 
homogeneous and isotropic to a sufficient level and the observer does not occupy a special 
location. (The small size of structures relative to the observed volume can be considered 
part of homogeneity and isotropy.) However, most studies of light propagation have 
neglected the influence of structures on the average expansion rate and the spatial 
curvature, assuming that they evolve according to the FRW equations. So the results 
should not be interpreted as evidence that the effect of structures is small, only that it 
is plausible that their effect can be considered in terms of the overall geometry. 

Even if that is true, one still needs to derive the proper distance as a function of 
redshift in terms of the expansion rate and spatial curvature, which describe the overall 
geometry. In general, the average spatial curvature does not evolve like a~^, where a is 
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the scale factor, and simply inserting spatial curvature with arbitrary time-dependence 
into the FRW relations does not make sense [68]. (Light propagation in terms of a 
scale factor which takes into account the influence of structures, but with the spatial 
curvature evolving like in the FRW case, was discussed in [69].) We will only study the 
evolution of the overall geometry as described by the scale factor and average spatial 
curvature, and will not consider the propagation of light. 

From the practical point of view, treatment of the passage of light in terms of a 
scale factor (even without any spatial curvature) has been very successful in accounting 
for the observations. If the observed deviation from the Einstein-de Sitter model was 
due to effects which cannot be expressed in terms of the average geometry, one would 
not expect different observations to agree so well with the simple scale factor treatment. 
For example, in the 'Hubble bubble' models, the density gradient which accounts for the 
luminosity distances to SNe la typically does not explain the angular diameter distance 
to the last scattering surface, or the scale of the baryon acoustic oscillations [47,67]. 

2.2. The expansion rate in a clumpy spacetime 

Inadequacy of the FRW description. Assuming that the influence of dumpiness on the 
passage of light can be encapsulated in the scale factor and the spatial curvature scalar, 
we are left with the question of how non-linear structures affect the expansion rate and 
the spatial curvature. 

The linearly perturbed FRW equations do not describe the local evolution of the real 
universe. Considering the linearly perturbed Einstein-de Sitter metric (and neglecting 
the decaying mode), the local expansion rate measured by a comoving observer is 
6 = 3H — 6H, where H is the background Hubble parameter in terms of the proper time 
and 5 cx a is the density contrast [25] . When 6 becomes of order unity, these relations fail 
to describe the real behaviour. In the case of underdensities, this is obvious, since the 
density contrast cannot decrease below —1. For typical spherical overdense regions, the 
difference between the linear and non-linear evolution is well-known from the spherical 
collapse model [70] (see [71-73] for reviews). 

It has been argued that even if there are large deviations locally, the statistical 
homogeneity and isotropy of the universe implies that the average evolution follows 
the FRW equations. However, a space which is statistically homogeneous and isotropic 
does not in general evolve on average like a space which is exactly homogeneous and 
isotropic. Exact homogeneity and isotropy leads to the FRW equations. Statistical 
homogeneity and isotropy simply supports the assumption that one can make sense of 
the observations in terms of the overall geometry. This does not imply that the scale 
factor follows the FRW equations, essentially because the spatial curvature in a clumpy 
spacetime in general evolves differently from the FRW case. (A test of whether the 
metric has the FRW form was proposed in [74], based on the specific FRW evolution of 
spatial curvature.) 

The influence of inhomogeneity and/or anisotropy on the expansion rate is known 
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as backreaction, and has been discussed in a number of papers [24,25,28,29,75-81]; 
see [29,62,82] for reviews and more references. (See [83-87] for discussion of backreaction 
during inflation.) 

Statistical homogeneity and isotropy. The early universe was highly homogeneous and 
isotropic, so the FRW equations provide a good approximation of the average evolution, 
and inhomogeneities and anisotropies are described by linear perturbations around the 
FRW model. Each hypersurface of constant proper time practically corresponds to a 
single expansion rate (with tiny variations). As perturbations grow and structures form, 
a given moment of time no longer corresponds to a single expansion rate. Instead, the 
hypersurface of constant proper time contains regions in different stages of expansion, 
some of them collapsing or static. On sufficiently large scales, the spatial distribution 
of expansion rates at each moment is statistically homogeneous and isotropic, and we 
can meaningfully call the local expansion rate averaged over all regions the expansion 
rate at that time. 

In order for the average expansion rate at a given time to be a useful quantity for 
describing the passage of light through structures, it is important that the expansion 
rate changes slowly compared to the time it takes for light to cross typical structures. 
Light has to have time to pass through different regions and sample the distribution of 
expansion rates before it changes appreciably. Otherwise, it does not make sense to use 
an expansion rate averaged over many regions: one would have to describe the passage 
of light through the individual regions. 

This condition is well satisfied in the real universe. For typical supersymmetric 
weakly interacting dark matter, the first structures which form around a redshift of 
2 ~ 40 — 60 have sizes of the order 10~^H^^ [88], and typical largest structures today 
have sizes around 10 /i~^Mpc ^ 3 x 10~^H~^ (where h parametrises the present-day 
Hubble rate, Hq = lOO/i km/s/Mpc). A more important quantity is the homogeneity 
scale, by which we mean the scale where averages converge to their asymptotic value. 
The fractal dimension of the set of galaxies around us indicates a homogeneity scale 
of 70-100 h~^Mpc [89,90], while studies of morphology suggest that it is at least 200 
h~^Mpc [91]\. In either case, the homogeneity scale is ~ 10~^H~^, so light rays pass 
through several representative regions of the universe in one Hubble time, which is the 
timescale for significant change in the expansion rate. 

A related concern about the applicability of the averaged expansion rate is the 
question of over which scale the averaging is done. To get a representative sample, 
the averaging scale should be at least as large as the homogeneity scale. Because of 
statistical homogeneity and isotropy, it should not matter if the averaging scale is larger 
than this. In practical terms, most cosmological observations of the expansion rate 
probe distances larger than the homogeneity scale, apart from SNe la at small redshifts. 

I The fact that there seems to be a significant contribution to our motion due to the Shapley 
Supercluster at a distance of 130-180 h~^Mpc [92] suggests that the homogeneity scale is on the 
larger side. Alternatively, our location in the universe may be rather untypical. 



Evaluating backreaction with the peak model of structure formation 



9 



For a statistically homogeneous and isotropic distribution of structures, varying the 
averaging scale as discussed in [93] should be relevant only for local observations, where 
the average scale factor treatment may be anyway problematic. (The observed large 
angle anomalies in the CMB might be related to local departures from the scale factor 
description; see [29] for discussion.) 

The issue of statistical homogeneity and isotropy is related to the choice of the 
hypersurface of averaging. The scale factor describes the volume of the hypersurface 
of proper time, but why take that hypersurface? Observations are organised on 
hypersurfaces of constant redshift, not proper time. However, if we can identify the 
redshift with a scale factor in the usual way, 1 + z = a{t)~^, then we have a one-to- 
one correspondence between the redshift z and proper time t, and the hypersurfaces of 
constant proper time and constant redshift agree. 

The question of the choice of hypersurface is also present in FRW models, and 
the hypersurface of constant proper time is selected because it is the hypersurface of 
homogeneity and isotropy. In realistic models, there is no such exact symmetry, but one 
still expects the hypersurface of constant proper time to agree with the hypersurface of 
statistical homogeneity and isotropy, and the argument for choosing this hypersurface is 
the same as in the FRW case. The evolution of structures proceeds according to proper 
time, so were one to tilt the hypersurface, different parts would contain structures 
which have evolved for different amounts of time, breaking statistical homogeneity and 
isotropy. 

The notion of statistical homogeneity and isotropy with a given homogeneity scale 
is easy to grasp intuitively: the universe consists of statistically identical boxes of a 
size given by the homogeneity scale. When evaluating any average physical quantities 
inside one box, the result should be independent of the location and orientation of the 
box, up to statistical fluctuations. This does not imply that the average quantities must 
be those of a model with exact homogeneity and isotropy. In terms of the Buchert 
equations discussed below, the backreaction variable Q saturates at the homogeneity 
scale, but not necessarily to zero. This view neglects long-range correlations, and in 
practice the degree of statistical homogeneity and isotropy increases as the size of the 
boxes grows, up to the level 10~^ given by the primordial perturbations. The concept 
of statistical homogeneity and isotropy in a spacetime with non-linear structures should 
be made more rigorous. Studies of the choice of hypersurface in backreaction [86,87,94] 
have not touched on this issue. 

3. The Buchert formalism 

3.1. The local equations 

The dust assumption. We assume that the matter content of the universe can be 
described as dust, i.e. an ideal fluid with zero pressure. This is not true on small scales 
(for example, the matter in the solar system cannot be treated as a pressureless ideal 
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fluid), so there is an implicit averaging involved, distinct from the large-scale averaging 
we are going to discuss. Essentially, the assumption is that on sufficiently large scales 
we can model a complex system of discrete small-scale structures as a continuum of 
infinitely fine particles. The validity of this approximation has been studied in the work 
on discreteness and the fluid approximation in N-body simulations [95]; indirect support 
can be found in the 'renormalisability' of Newtonian gravity in [96] (see also [97]). 

The Einstein equation. For dust, the Einstein equation reads (taking the cosmological 
constant to be zero) 

Gal3 = 87rG-!<jTaf^ 

= SuG^^pUaUp , (1) 

where Gap is the Einstein tensor, Gn is Newton's constant, Taf3 is the energy-momentum 
tensor, p is the energy density and is the velocity of observers comoving with the 
dust. 

Following the covariant coordinate-independent approach, the Einstein equation (1) 
can be decomposed into scalar, vector and tensor parts. This is a local decomposition 
with respect to general coordinate transformations, not a global decomposition in terms 
of the symmetry of a background, as in perturbed FRW spacetimes. (We are not 
assuming any symmetries, or making a division into background and perturbations.) 
For reviews of the covariant approach, see [98-101]. 

We want to discuss averages. Since vector and tensor quantities cannot be 
straightforwardly averaged§, we will consider the scalar part of the Einstein equation. 
The Einstein equation (1) has two scalar components, and the covariant conservation 
law yields a third equation [98,99,104,105] (see e.g. [77] for the full system of equations): 

9 + -e^ = - AttG^p - + 2c^2 (2) 
3 

V =87rG^p-^f^R + a'-u;' (3) 

P + Op =0, (4) 

where a dot stands for derivative with respect to proper time t measured by observers 
comoving with the dust, 6 is the expansion rate of the local volume element, cr^ = 
|(T"^(Ja/3 > is the scalar built from the shear tensor aa/3, oj^ = \uj"'^ujai3 > is the 
scalar built from the vorticity tensor uJap^ and '-^^i? is the Ricci scalar on the tangent space 
orthogonal to the fluid flow. The acceleration equation (2) is known as the Raychaudhuri 
equation, and (3) is the Hamiltonian constraint. The equations are exact, and valid for 
arbitrary large variations in density, expansion rate and other physical quantities. 

Vorticity. When the vorticity is zero, the tangent spaces orthogonal to the fluid flow 
form spatial hypersurfaces which provide a foliation that fills the spacetime exactly once. 

§ Though see the work on the Ricci flow [79] and the 'macroscopic gravity' formahsm [102]. On the 
relation of the latter to the averaging scheme used here, see [103]. 
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These flow-orthogonal hypersurfaces are also hypersurfaces of constant proper time of 
comoving observers. If the vorticity is non-zero, no such hypersurfaces exist. (See [105] 
for discussion of the Ricci scalar of the tangent spaces in the case of non-zero vorticity.) 

The covariant approach deals directly with physical quantities, and there is no need 
to introduce coordinates. Nevertheless, it can be shown by explicit construction that if 
we wish to choose coordinates, it is always possible to take goo = —1, and simultaneously 
obtaining goi = is possible if and only if the vorticity is zero (here stands for time 
and i for the spatial directions) [106]. For irrotational dust, the synchronous metric can 
be adopted locally without any loss of generality. In the covariant formulation, there 
are no artifacts related to the choice of coordinates, since there are no coordinates, 
and all variables are physical observables. However, there would be no problem (when 
the vorticity is zero) in using the synchronous comoving coordinates. The results for 
physical quantities of course do not depend on the coordinate system. (For example, in 
perturbative backreaction calculations one obtains the same results in the longitudinal 
and comoving synchronous gauges [81].) 

We take the vorticity to be zero. For dust, if vorticity is initially zero, it will remain 
zero. Furthermore, vorticity in linearly perturbed FRW models corresponds to vector 
perturbations, and it decays with expansion (unlike shear). However, the description 
of matter as a pressureless ideal fluid will break down on small scales when non-linear 
structures form, due to shell-crossing. After any vorticity is generated in gravitational 
collapse [107], it will be amplified by the collapse, and will formally diverge at the same 
time as the density [108]. A static dust structure, with ^ = 0, is possible only if large 
amounts of vorticity are present, as (2) shows: vorticity has to balance the contribution 
of both the energy density and the shear on the right-hand side. For stabilised dust 
structures, vorticity is the dominant contribution. (As the approximation of treating 
matter as dust breaks down, effects such as velocity dispersion and pressure can be also 
important in the stabilisation [109].) 

Nevertheless, just as we assume that the small-scale breakdown of the description of 
matter as dust is not important, we assume that the vorticity inevitably present at small 
scales can be neglected in discussing the overall large-scale evolution. (This assumption 
is also involved in the usual perturbed FRW treatment, with rotationless ideal fluids and 
a well-defined cosmic time.) Vorticity is probably only relevant in stabilising structures. 
Since we do not need to consider the details of stabilisation in our calculation, we do 
not expect the assumption of zero vorticity to be important. 



3.2. The average equations 

Defining the average. We follow the formalism introduced by Buchert in [76,77,80]. 
The spatial average of a quantity is its integral over the hypersurface of constant proper 
time t, divided by the volume of the hypersurface 

{fm^^, (5) 
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where Cq/^^ = rjap-fsu^ is the volume element on the hypersurface of proper time, r]ai3'y& 
being the spacetime volume element. 

The scale factor is defined simply as the volume of the hypersurface of constant 
proper time to power 1/3 

(a) 

where a has been normalised to unity at time to? which we take to be today. As 9 is the 
volume expansion rate, this definition of a is equivalent to 3a/a = {6). We will also use 
the notation H = a/a. 

The Buchert equations. Let us take the average of the equations (2)-(4). The resulting 
Buchert equations are [77]: 

3- = -47rGN(p) + Q (7) 
a 

3^ = 8vrGN(p)-^((^)i?)-iQ (8) 

dt{p) + 3-{p) = 0, (9) 
a 

where the backreaction variable Q contains the effect of inhomogeneity and anisotropy: 

Q^lm-{e)')-2{a') . (10) 

The integrability condition for the average Raychaudhuri equation (7) and the average 
Hamiltonian constraint (8) is 



d,e^R)+2-{^'^R) = -d,Q-6-Q, (11) 



a a 
The Buchert equations (7)-(9) are exact for the averages when matter consists of 
irrotational dust. (The corresponding equations in the Newtonian case were derived 
in [76], and the case with non-zero pressure was considered in [80].) The variance 
of the expansion rate in Q is a new term compared to the local equations. It has 
no counterpart in the local dynamics, and may be called emergent in the sense that 
it is purely a property of the averaged system. If the variance is large enough, the 
average expansion rate can accelerate, even though the local expansion rate decelerates 
everywhere. (In physical terms, the average expansion rate can grow because the volume 
occupied by faster expanding regions rises.) 

The density parameters. As in the case of FRW models, we can parametrise the 
different contributions to the expansion rate with relative densities. Dividing (7) and 
(8) by SH"^, we have [77,110] 

q= -^- = ln^ + 2nQ (12) 
H'^ a 2 

i = n^ + nR + nQ , (13) 
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where Q^a = SnGNip) /{SH'^), Qr = -{'^^^R)/{6H^) and Qq = -Q/{6H^) are 



the density parameters of matter, spatial curvature and the backreaction variable, 
respectively. As seen from the definition of Q in (10), the backreaction density parameter 
is just minus the relative variance of the expansion rate, plus the contribution of shear: 



The backreaction density can have either sign, and can become arbitrarily 
negative. Correspondingly, q can cross the value —1, and the contribution Qr + Qq can 
change sign. These features are not captured in the parametrisation of backreaction in 
terms of a scalar field, the 'morphon' [111]. 

3. 3. First integrals of the Buchert equations 

First integral in terms of Q. The scalar parts of the Einstein equation (2) and (3) 
together with the conservation law (4) do not form a closed system on their own. (In 
particular, the propagation equations for the shear and vorticity tensors cannot be 
reduced to scalars.) There are four unknowns (9, p, a"^ — cu^, ^^^R) and three equations. 
The Buchert equations (7)-(9) are similarly underdetermined, with four unknowns 
(a, (p), Q, i^^^R)) and three equations. 

If we give as extra input the evolution of one of the variables (or a relation between 
the variables), the average system is completely determined. In particular, we can find 
the average expansion rate from the average spatial curvature {^^'^R) or the combination 
of variance and shear given in Q. We can make this explicit with the first integral of 
the Buchert equations (7)-(9). 

Taking into account that the conservation of mass (9) implies (p) oc a~^, we can 
integrate (7) to obtain 



where K is an integration constant related to the spatial curvature and (po) is the 
average energy density at time to- When backreaction vanishes, we recover the FRW 
Hubble relation, and the spatial curvature is 6Ka~^, as usual. When the spacetime 
is clumpy and backreaction is present, Q ^ 0, the average spatial curvature evolves 
non-trivially: 




(14) 




(15) 



First integral in terms of {^^^R) . It is also instructive to write the first integral in terms 
of the spatial curvature. Expressing things in terms of {^^^R), we have from (7)-(9), 



J a' ^ ' ^ ^ 



where C is an integration constant. The backreaction variable is 




(17) 
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When the average spatial curvature evolves like a~^, the C term gives the 
backreaction variable Q, which in accordance with the integrability constraint (11) then 
evolves like a~^. We can make this more explicit by decomposing the average spatial 
curvature as {^^^R) = 6Ka~^ + A{^^^R). The backreaction variable is|| 



In general, underdense regions are negatively curved and expand faster than 
the average, while overdense regions are positively curved and expand slower. The 
evolving distribution of these regions determines the average spatial curvature. At 
late times, one would expect the faster expanding regions to dominate the volume, 
and the spatial curvature to be negative. The intertwining of the expansion rate and 
the spatial curvature is quantified in (19). While Q is non-local, {^^^R) is simply the 
average of a local quantity; they are interchangeable via the integrability condition 
(11). Understanding backreaction in terms of an average over the local spatial curvature 
will be useful when we discuss the importance of non-Newtonian aspects of gravity in 
backreaction in section 5.2. 

We need only know the evolution of the variance of the expansion rate minus the 
shear in order to reconstruct the complete expansion rate (up to a constant related to 
the spatial curvature). There is no need to specify the full metric, only some statistics. 
Indeed, it is unfeasible to write down a metric that would be a realistic description of 
the structures in the universe. Exact solutions such as the spherically symmetric LTB 
model are useful for illuminating specific aspects of backreaction such as the choice of 
hypersurface [94] and demonstrating acceleration unambiguously [30,113,114]. However, 
exact solutions, even ones with no symmetry such as the Szekeres model [115], are 
limited in scope. While more complicated than the LTB model, they are not much 
closer to a realistic picture of the universe with its hierarchical layers of individually 
complex but statistically homogeneous and isotropic structures. In order to evaluate 
the importance of backreaction in the real universe, we need statistical knowledge about 
complex configurations of dust, not exact information about simplified models. 

II In [112] it was argued on the basis of a 2+1-dimcnsional study that the effect of backreaction is small. 
However, since the Ricci scalar in two dimensions is a topological invariant, the average two-dimensional 
spatial curvature is necessarily proportional to a~^, and the evolution of Q is trivial. Therefore the 
situation in 2+1 dimensions is not representative of the 3+1-dimcnsional case, where spatial curvature 
is dynamical. There were also two other arguments presented in [112], one based on second order 
perturbation theory and the other on a solution with parallel walls. The first calculation is incorrect, 
because it assumes that the background energy density is the same as the average energy density, which 
is not true beyond the linear level. In the second case, variations in the local expansion rate are small, 
so it is irrelevant for the real universe. 





and the Hubble equation (16) reads 




(19) 
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4. Backreaction with the peak model 
4.I. Setting up the model 

The universe divided into regions. We want to have a model of the universe that 
does not violate large-scale statistical homogeneity and isotropy, and includes a realistic 
ensemble of evolving structures. We will treat each structure as an isolated region, and 
take the nested nature of hierarchical structure formation into account in the evolution 
of the number density of the regions. 

We divide the hypersurface of constant proper time into non-overlapping regions 
labelled by (5. The fraction of proper volume occupied by each region is ^^(t) = /^^e/ f^e, 
where ^ is the integral over region 6 at time t. This division is completely general. We 
now take the regions to be isolated and spherically symmetric. This is an extension of 
the two- region toy model studied in [28,29] to cover a realistic distribution of structures. 
An ensemble of only spherically symmetric regions cannot fully cover the hypersurface of 
proper time (and does not form a connected space), and the model has to be understood 
in a statistical sense. (We come back to this point in section 5.2.) 

We treat the spherical regions with Newtonian gravity. Their average evolution is 
then given by the spherical collapse model (and its underdense equivalent), according 
to which they evolve like the corresponding FRW universe [70] (see [71-73] for reviews). 
In terms of the Buchert equations, this follows from the result that in Newtonian 
gravity, Q vanishes for spherical symmetry, so the Buchert equations reduce to the 
FRW equations [26,78]. Expansion of a region with positive density contrast slows down 
more as the perturbation grows, until it turns around and collapses, finally stabilising 
at a finite size and density. Inside a region with negative density contrast, expansion 
decelerates less as the region becomes emptier. 

We take the individual regions of the ensemble to represent structures at a definite 
state of expansion or collapse. In other words, the label 6 corresponds to a value of 
the average expansion rate in a region. In the spherical collapse model, the expansion 
rate is in one-to-one correspondence with the linear density contrast (hence the letter 
6), which is the density contrast that the region would have relative to an Einstein-de 
Sitter universe if its evolution had continued as in the linear regime. 

The average expansion rate is 



The fraction of proper volume in structures with linear density contrast 6 at time t 
has been decomposed into two parts, vs{t) = ssf{S,t)/{J^^d6ssf{S,t)). (Note that 
the linear density contrast can be arbitrarily negative, unlike the real density contrast, 
which is bounded from below by —1.) 

The first term ss = asit)^ / aEdsit)^ is the volume of a region with linear density 




JZ,d6ssf{6,t)Hs{t) 
JZdSssf{6,t) 



(20) 
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contrast 6 relative to the volume of the Einstein-de Sitter universe, where OEds ^ ^^"^ (we 
will also use the notation -ffEds = 2/(3t)). It is due to the difference in the expansion 
rate between the different regions: more underdense regions expand faster and have 
therefore grown larger. (The scale factor OEds has been introduced out of convenience. 
It appears both in the numerator and the denominator and does not depend on 5, so it 
does not affect the average.) 

The second term f{S,t) is the fraction of the initial volume in regions with linear 
density contrast 6. This is updated to the present volume by the first term. If we 
completely ignored interactions between regions, f{S,t) would be directly given by the 
linear spectrum of perturbations. However, we want to take into account the nested 
nature of cosmological perturbations, and the merger of structures into larger entities. 
As overdense regions collapse and stabilise, they form larger structures which in turn 
slow down and collapse, and underdense regions have similar hierarchical evolution. We 
will include this feature using the peak model of structure formation. 

The premise of the peak model of structure formation is that structures are 
identified with maxima of the linear, Gaussian density field, smoothed on some scale 
R [116]. In the original application, all peaks above a fixed density threshold were 
considered to be stabilised non-linear structures. We will use the peak number density as 
a function of density contrast as the number density of isolated regions having that linear 
density contrast. (Since the density field is Gaussian, the distribution of underdense 
troughs is the same as the distribution of overdense peaks.) The number density can be 
converted into the fraction of mass. Since the early universe was very smooth, this is 
also the fraction of the initial volume which has ended up in structures of a given linear 
density contrast, in other words our f{S,t). 



The backreaction variable. Our treatment will give us the average expansion rate 
directly, without needing to go via the backreaction variable Q. However, interpreting 
the expansion rate in terms of the Buchert equations (7)-(9) and the density parameters 
(12) will help to understand its evolution. The backreaction variable for the model is 



2 
3 

2 
3 



^2^ 

32\ ' ' 



oo \J — oo 



d5vs{0')s- / d5vs{e)s -2 / d6vs{a')s 



oo / POO \ 2 

d6vs{e)^s- / d6vs{9)s 



oo \J — oo 

oo POO 



+ I ( / d6vs{e^)s - I d6vs{e)l ) - 2 / d6vs{a^)s 

'•^ \J —oo J —oo / J —oo 

POO / POO \ 2\ noo 

/ d6vsH]-( d6vsHs) +/ d6vsQs, (21) 

J —oo \J —oo J I J —oo 

where is the average over region 6. The total backreaction variable Q is not the 
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sum of the regional backreaction variables, due to the non-local variance term. Since 
we consider spherical regions using Newtonian gravity, we have Qs = [26,78], and 
the overall backreaction variable can be calculated from the regional average expansion 
rates. 



The spherical collapse/expansion model. The expansion rates of the regions follow the 
spherical collapse model (and its underdense equivalent), which we summarise here 
(see [73] for another useful summary). 
For overdense regions we have 



sm0(0 — sm< 

Hx+t - 



(1 — cos (f)y 
_ al+ 2(1 — cos I 
~ (^Eds 9(0- sin < 

1 — COS 



5+ = A62/^(0-sin0)^/^ (22) 

where (p is the development angle which runs from to 27r and is the linear density 
contrast. The structure collapses to a singularity at 27r, so we take the expansion rate 
to stabilise discontinuously to zero at = 2-7r. We set the volume, expansion rate and 
spatial curvature of regions with a linear density contrast larger than 6~^{27i) ^ 1.7 to 
zero. It would be more realistic to set the volume to a finite constant and the spatial 
curvature to some constant proportional to the stabilised energy density, but this makes 
little difference for our purposes. In our calculation, each region is weighted by its 
relative volume, so the contribution of collapsing regions goes to zero at the collapse, and 
the details of stabilisation do not matter. (For example, the results would be unchanged 
if we stopped the evolution at ^ = 37r/2 instead of 27r, as is sometimes done.) While 
the expansion rate is discontinuous at the stabilisation, the volume-weighted expansion 
rate ss+Hg+t is continuous. The combinations Hs+t, ss+ and {^^^R)s+t'^ do not depend 
explicitly on t, they are functions of alone. 

Correspondingly, for underdense regions we have 

sinh 0(sinh (f) — (j)) 



H,-t 



(cosh0 — ly 
a^- _ 2(cosh0- 1)3 



a-'^Ac 9(sinh — 0)^ 

A\2 



EdS 

(sinh ( 



{^'^R)s-t' = -6 



(cosh — ly 



5- = -— e'/^^sinh - 0)2/3 , (23) 
20 

where the development angle runs from to oo, corresponding to increasing time and 
decreasing density contrast. 
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Regions with zero density contrast have H-^dst = 2/3, s^o = 1 and zero spatial 
curvature. 

Peak statistics. For a hnear density field with Gaussian, statistically homogeneous and 
isotropic perturbations smoothed on scale R, the number density of peaks (troughs) of 
height (depth) u = 6/aoit,R) is [116] 



n{u, R) = e 2' 



(27r)2i?,(/?)3 
where the function G(z/, 7(_R)) is 



G{u,^iR)) 



I -(^-7kl)^ 

/o V27r(l - 72) 



(24) 



(25) 



with 



F{x) 



X 



3x 



erf 



+ erf - - 



(svr) 



31x2 8 



5'^ 




(26) 



The functions 7(-R) and R^{R) are defined as (note that the explicit time-dependence 
in aj(t,R) cancels out) 

af{t,R) 



7(i?) 



^a'o{t,R)ai{t,R) 
'ai{t,R) ' 

where the spectral moments cr'^it, R) are 



-R*(-R) 



(27) 



a]it,R) 



dk 
T 



4 



k^Aj{k,t)T{kYW{kRf 



dk 



9 (ctEdS-f^Eds)^ 

4 

9 (aEds-f^Eds)^ 



^k'+^Al{k)T{kfW{kR)^ 
" ^A;^+4T(A;)2iy(A;i?)2 , 



f28) 



where is the primordial power spectrum of density perturbations, is the primordial 
power spectrum of metric perturbations, taken to be scale-invariant with amplitude 
= {3x 10~^)2, T{k) is the transfer function and W{kR) is the window function, taken 



to be Gaussian, W{kR) 



(When going from metric perturbations to density 



perturbations, we have neglected the long-wavelength modes, as we are interested in 
sub- horizon perturbations.) The linear density field is taken to evolve in an Einstein-de 
Sitter universe. We are neglecting the 'backreaction of backreaction', the effect of the 
change of the average evolution on the perturbations. One would need to take this into 
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account to see if there are, for example, oscillations in the expansion rate as suggested 
in [28,29]. For work on perturbations in a backreaction context, see [27,29,78,117-119]. 

To go from the peak number density to the fraction of mass (or initial volume), we 
need to specify the mass associated with each peak. We multiply the number density 
(24) by the volume under the Gaussian window function, (27r)^/^i?^, and introduce a 
constant to account for the fact that the total mass is not correctly normalised. The 
fraction of mass in peaks of height u = 6/aQ{t, R) is then 

The total mass defined this way is not constant with R. We fix the normalisation 
constant by demanding that asymptotically all mass resides in peaks or troughs, 

dz//(z/, i?) — i> 1 as i? — i> oo. This gives ~ 1.97. In this treatment, all peaks 
(and troughs) contain the same amount of mass, regardless of their height (or depth). 
A normalisation factor which depends on u would probably be more appropriate, but 
we want to keep the treatment simple. The fraction of volume which is not in peaks 
or troughs is taken to expand like the Einstein-de Sitter universe, H^^st = 2/3. For 
discussion of mass assignment for peaks and troughs, see [73,120,121]. 

One factor which we have not taken into account is the peak-in-a-peak problem. 
The distribution function of peaks (29) does not take into account that lower peaks may 
be submerged in higher ones. There is an equivalent trough-in-a-trough problem, and 
perhaps most importantly, the trough-in-a-peak problem, as some underdense regions 
are extinguished by larger overdense regions. There is no corresponding peak-in-a-trough 
problem, and this asymmetry transfers mass from the underdense to the overdense 
regions [73]. Our treatment does not include this effect, and we have equal mass in the 
underdense and overdense regions, since the density field is Gaussian. 



Smoothing and time evolution. We determine the smoothing length R by fixing crQ(t, R) 
to a given value at all times. We take this value to be unity, (TQ{t, i?) = 1, so we have 
u = 6. Since the linear density contrast evolves in time, crQ{t, R) oc a^^g, the smoothing 
length R{t) will also evolve, growing with time. This in turn translates into evolution of 
the distribution function f{S,R{t)). The smoothing scale R{t) can be regarded as the 
typical size of structures which are forming at time t. 

All time evolution is taken into account in this statistical manner in the distribution 
function f{6,t), since the individual regions are by definition fixed at a given state 
of expansion or collapse. This averaging of regional values of the expansion rate is 
somewhat different from the Buchert formalism, where the basic quantity is the local 
volume element. While the density parameters defined in (12) can be used to understand 
the state of the universe at each moment, the evolution of the expansion rate is not 
completely captured by the Buchert equations, since the indirect treatment of the local 
evolution of mergers by smoothing goes outside the dust approximation. 

In physical terms, smoothing with a window function involves the assumption that 
we need not consider the details of structures at the smoothing scale. It is a simplification 
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Figure 1. The BBKS (blue, solid) and BD (red, dasiied) transfer functions as a 
function of fc/fccq. 



for the merger of structures, both overdense and underdense, into larger entities. In the 
present setting, this is part of our assumption of replacing continuous evolving structures 
by disjoint spherical regions. In a full description of structures, there would be a well- 
defined density contrast and expansion rate at each point. When we model the universe 
with regions which are each associated with a regional average density contrast and 
expansion rate, smoothing is introduced. 



The transfer function. The statistics of structures in the peak model are determined 
by the functions 7 and R^, defined in (27). With a fixed primordial power spectrum, the 
behaviour of these functions is given by the transfer function. 

We will consider two approximations for the CDM transfer function. The BBKS 
transfer function [116] is 

m = H^+2M,) 

2Mq [1 + 3.89g + (16.1g)2 + (5.46g)3 + (6.71g)4]^/^ 

with q = fce^V(13.7fceq) ( assuming three massless neutrino species), where /b = 
^b/^m is the baryon fraction [122] and ~ 13.7ci;~^ Mpc is the wavelength of 
the perturbations that enter the horizon during the matter-radiation equality. Here 
= ^moh"^- The BBKS transfer function is a fitting formula to numerical results, and 
the large k limit has been also analytically derived [122]. We fix the baryon fraction 
at /b = 0.2. The results do not significantly depend on /b directly, but the BBKS 
approximation of the numerical transfer function calculated with CAMB becomes worse 
with increasing /{,; for /b = 0.2, the error at large k is 20-30%. 

For comparison, we will use the simple transfer function introduced in [41] by 
Bonvin and Durrer, which we will call the BD transfer function, 

= . , on n ^. ' (31) 



P{k/ kcq^ 



where /5 = 3 x 10"! 

The BD transfer function is a simple approximation to the behaviour of the 
evolution of CDM perturbations captured in more detail by (30). It is not quantitatively 
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accurate, but it will help to estimate the dependence of the results on the transfer 
function. 

Qualitatively, the BD transfer function has the same behaviour as the BBKS 
transfer function. Perturbations with wavenumber much larger than fceq enter the 
horizon deep in the radiation dominated era, so they are damped by while 
perturbations with much smaller wavenumber enter in the matter-dominated era, so 
their amplitude is unsuppressed. The evolution in the slope of the BBKS transfer 
function is shown in figure 1. Both the BBKS and the BD transfer function are missing 
the cut-off at small scales due to free-streaming [88, 123]. However, we will see that the 
small-scale behaviour is not important for our results. 

4-2. Preferred time in structure formation 

The size of structures. Before evaluating the average expansion rate using the peak 
model, we will briefly discuss the preferred time that is involved in the formation of 
CDM structures from a scale-invariant spectrum of primordial perturbations. We would 
expect to see this timescale reflected in the evolution of the expansion rate, regardless 
of the statistical model used to describe structures. 

Because of the change of the slope of the CDM transfer function, the size of 
structures relative to the Hubble size saturates around 10-100 billion years [28,29]. 
This can be seen as follows. Typical structures forming at time t have size R{t) defined 
by cr^it, R) = 1. Inverting (28) to solve for R and dividing by the Hubble radius, we get 
the relative size R/{aH)~^ as a function of time. 

The timescale is determined by the time of matter-radiation equality teq ~ 1000c<j~^ 
years. This value assumes three massless neutrino species with abundances determined 
by thermal equilibrium, but it is independent of late-time cosmology, as long as the 
energy density of matter evolves like and the energy density of radiation evolves like 
a~^. In a clumpy space, the first follows from the conservation of mass (9). There is no 
such conserved quantity for radiation, so in general the energy density of radiation does 
not necessarily evolve like [80]. However, if the number density of radiation quanta 
is conserved and their change in energy (i.e. redshift) is, at least on average, related to 
the scale factor hj 1 + z = a~^, the radiation energy density will be proportional a~^. 
This relates the studies of the passage of light to the study of the Buchert equations in 
the case of non-zero pressure [80]. 

We adopt the value uj^a = 0.1 for all plots, giving tgq ~ 10^ years. This value of u^a 
could be reasonably moved up or down by a factor of 2, so the timescale in the plots 
could be shifted by a factor of 4 in either direction. Model-independent estimates of 
Qrao can be summarised as 0.15 > Q^ao ^ 0.35 [124]. Regarding the Hubble parameter, 
the value from SNe la observed with the Hubble Space Telescope has been quoted as 
h = 0.72 ± 0.08 [125] or h = 0.62 ± 0.05 [126], depending on the treatment of Cepheids 
(see section 5.1). For comparison, the model-dependent value determined from fitting 
the ACDM model to the WMAP3 data is = 0.127toml P]- 
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(a) (b) 

Figure 2. The size of structures R/{aH)^^ relative to the asymptotic size as a function 
of time, in biUions of years, for (a) the BBKS transfer function and (b) the BD transfer 
function. 



At early times, the damping due to radiation domination cancels the 

k"^ enhancement of the density perturbations relative to the scale-invariant metric 
perturbations. Therefore the spectrum of density perturbations depends only weakly on 
the scale, so structures on different scales collapse nearly at the same time. The initial 
structures are small and grow rapidly. Beyond keq, the amplitude is not suppressed, so 
after the wavenumber of collapsing perturbations reaches kg^, there is no scale in the 
system anymore, and the relative size of the structures is constant. The asymptotic size 
is roughly \/A ^ 5 x 10^^. 

We show the value of R/{aH)~^ relative to its asymptotic value in figure 2 for the 
BBKS and BD transfer functions. The timescale of the evolution of the size is sensitive 
to the behaviour of the transfer function at small wavenumbers. For the BBKS transfer 
function, the relative size enters the saturation regime at some tens of billions of years. 
For the BD transfer function, this happens at a few billion years or so. In both cases, the 
evolution in the size practically saturates well before the perturbations with wavenumber 
kcq become non-linear, which happens around 3500 billion years for the BBKS transfer 
function and 2000 billion years for the BD transfer function. 

The rough agreement of the era when structures reach their maximum size with 
the time when acceleration has been observed (about 10 billion years) is encouraging 
from the point of view of the coincidence problem. A rough way to understand the 
preferred time in terms of the time of matter-radiation equality and the amplitude of the 
primordial perturbations is as follows. The era when the structures with wavenumber 
keq form is determined by the condition a^{t,k~^) = 1. Approximating the window 
function W{kR) with the step function, the transfer function with unity for k < k^q 
and recalling that kcq = {aH)~^, we get from (28) the time t = {S/Ay^Heq ~ 3000 
billion years. The approach to the asymptotic value of the transfer function is slow, 
and the saturation regime starts a couple of order of magnitude earlier than this. (It 
is apparent in the BD transfer function (31) that the turnover scale is ~ 10/Ceq rather 
than keq- Naively using 10/ceq in the previous argument would bring down the time to 
3 billion years, since t oc (ctEds-f^Eds)"^-) Exactly because the size changes slowly, it is 
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somewhat arbitrary what one calls the preferred time. Around 10-100 billion years the 
relative size of structures changes from growing rapidly to being almost constant, but 
one cannot be more precise than that. 

We will now get back to our model with the peak statistics, where the expansion 
rate is determined as a function of time without ambiguity, and see how it indeed changes 
around the saturation era. 



4-3. The expansion rate 

Increasing Ht. The expansion rate is given by 

/oo 
d5 vs{t)Hst 
-oo 

IL d^" ^s-f{5-, R)Hs-t + /,°° d^+ ss.fi5+, R)Hs.t + |(1 - f^^ d5 f{5, R)) 
/I ss^ f{6-, R) + d5+ R) + {1- IZ /(^' R)) 

where ss± and Hs±t are given in (22) and (23), and / is given in (29). The expansion rate 
is completely determined, there are no free parameters to adjust (unless one counts the 
baryon fraction in the BBKS transfer function). The average spatial curvature {^^^R) 
and the backreaction variable Q given in (21) are calculated the same way. 

We first show Ht as a function of the size of structures relative to the equality scale, 
r = keqR, in figure 3. Today (jQit, R) is unity on the scale of 8 h~^Mpc or slightly below, 
so given k~^ ^ 13.7u;~^ Mpc, we have r !^ ~ 0.05-0.2, placing the present time in 
the transition region. The transition era in the expansion rate is clear and the change 
looks rapid, particularly for the BD transfer function. However, because the growth 
of r as a function of time slows down as time goes on, the evolution is less steep as a 
function of t. 

In figure 4 we show Ht as a function of time. The behaviour is qualitatively similar 
for the BBKS and BD transfer functions. At early times, Ht is close to the FRW value 
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2/3, slightly higher because there are structures present. With the onset of saturation 
(i.e. the turnover in the transfer function), Ht rises to a value somewhat less than unity. 

The physical reason for this evolution is that at early times the volume occupied by 
the structures is small, so their impact is small. As the volume occupied by structures 
grows (along with the density contrast of typical structures), the expansion rate becomes 
dominated by voids, since their volume is large. This picture is in qualitative agreement 
with excursion set studies of voids [73, 127]. If all volume was in voids that were 
completely empty, we would have Ht = 1. Because the voids are not completely 
empty, and because there are overdense regions, the expansion rate asymptotes to a 
somewhat smaller value. This evolution could be expected on general grounds, as 
discussed in [29], but it is not trivial that the timescale comes out correctly, in agreement 
with the argument related to the size of structures presented in the previous section. 

As the Buchert equations (7)-(9) show, the effect of perturbations on the average 
expansion rate does not necessarily become large when they first become non-linear. 
The criteria for the breakdown of the perturbed FRW equations as a description of the 
local evolution and as a description of the average evolution are different. The latter 
breaks down only when non-linear density perturbations occupy a sizeable fraction of 
space^. 

No acceleration. In figure 5 we show the deceleration parameter q = —a/{aH'^) as a 
function of time. Even though Ht rises, there is no acceleration. The reason is that 
only the underdense regions are important, the overdense regions play almost no role. 
As a result, Ht does not slow down before rising. It is likely that in order to have 
acceleration, the overdense regions should first slow down the expansion rate, so that 
the relative difference is larger, as in the toy model discussed in [28,29]. In FRW models, 

% The time when backreaction becomes important was argued in [128] to be around today, by estimating 
when the perturbative corrections to the average expansion rate become of order one. However, the 
analysis does not take into account that the average decomposes into a product of two-point functions, 
and that spatial derivatives have to occur in pairs. The correct order of magnitude for the general term 
after eq. (8) in [128] is 10"^(a/ao)((5^)*"^^''/^. This does not select out the present day. 
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(a) (b) 

Figure 6. The density parameters (red, solid), (blue, dotted) and (green, 
daslied) as a function of time for (a) the BBKS transfer function and (b) the BD 
transfer function. 



there is a clear qualitative difference between acceleration and deceleration, related to 
the dominant energy condition. Here the situation is different: acceleration is just 
a quantitative question of the steepness of the Ht curve and the size of the density 
parameter Qq, and there is no principle involved. 

In figure 6 we plot the density parameters Q^, and Qq defined in (12). For 
both transfer functions, the backreaction density parameter Qq is small, below 0.04 in 
absolute value. (The backreaction contribution peaks in the transition era when Ht rises 
rapidly, but the change is too small to appreciate in the plots.) As noted in [78,129], the 
system can evolve far from the initial near-FRW behaviour, even though Qq is small at 
each moment. The backreaction density parameter Qq is a measure of the distance of 
the model from the FRW case. The smallness of Qq means that at each moment, the 
system is near some FRW model and evolves slowly between different near-FRW models. 
(Note that the asymptotic value of Qq is non-zero, so unlike an open FRW universe, 
the model does not become emptier without limit: structures are present, even when 
the volume is dominated by voids.) In order for the behaviour at some instant to be 
very different from a dust FRW model, \Qq\ has to be large. A backreaction variable 
\Qq\ > 0.2 is needed to obtain acceleration of the observed magnitude [29]. 

While the timescale of the change in the expansion rate comes out roughly in 
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agreement with observations, the effect of dumpiness is too small to give acceleration. 
Let us now discuss the assumptions involved in the model, and how to improve the 
treatment. 

4.4- Improving the model 

The spherical collapse/expansion model. Our model involves two parts: structures are 
treated as an ensemble of regions whose distribution is given by the peak model, and the 
evolution of the individual regions is taken to follow the spherical collapse/expansion 
model. 

Let us discuss the collapse/expansion model first. While the spherical collapse 
model works surprisingly well for the statistics of structures, it does not give an 
accurate description of the collapsing phase. In particular, there is no treatment 
of the stabilisation which ends the collapse. An effective treatment of the terms 
responsible for stabilisation was given in [130], and a improved model which covers the 
whole evolution from the linear regime to the stabilised phase was presented in [131]. 
One can also generalise into ellipsoidal collapse, and take shear and tidal effects into 
account [107,121,132]. The spherical collapse model in a backreaction context has been 
discussed in [119,133]. 

In general, a given value of the linear density contrast does not correspond to a 
fixed value of the expansion rate, unlike in the spherical collapse model. In terms of the 
Buchert equations (7)-(9), the backreaction variable Q is non-zero, and the equations 
do not have a unique solution. However, one could generalise the present treatment to 
cover more realistic structures by adding one layer of ensemble. One could integrate 
over a distribution of initial conditions to obtain the mean expansion rate corresponding 
to structures of a given density contrast, using the Buchert equations, as presented 
in [78, 133]. While the approximation that the expansion rate depends only on the 
density contrast (used in [130,131]) is probably not valid for a single structure, it is true 
by construction for the expansion rate averaged over different structures. 

There is no stabilisation problem with the underdense equivalent of the spherical 
collapse model, and the evolution smoothly approaches the empty case Hs-t = 1. 
However, the model breaks down if shells from inner regions of the void (which have a 
smaller density and therefore expand faster) catch up with outer shells, even though this 
is not apparent in the average description. When this happens depends on the initial 
density profile. For a top-hat density profile, shell-crossing occurs at linear density 
contrast S~ ~ —2.8, corresponding to H^-t = 8/9 [73]. In real structures, the shell- 
crossing singularity corresponds to the formation of a dense wall surrounding the void. 
After wall formation, it might be more realistic to take the expansion rate to be that of 
an outward-moving shell, Hs-t ^ 0.8 [134, 135]. Like in the case of overdense regions, 
the assumption of spherical symmetry may be questionable. While isolated voids grow 
more spherical [135,136], the shapes of real voids are affected by surrounding structures. 
Voids in simulations have complicated shapes [137,138], though the overall structure 
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resembles a foam of spherical voids [73]. The shapes of underdense structures could 
be taken statistically into account with an ensemble of initial conditions, as with the 
overdense regions. 

However, refinements to the collapse/expansion model are secondary to the issue 
of distribution of the overdense and underdense regions. For example, redoing the 
calculation above with the improved collapse model presented in [131] makes practically 
no difference to the results. The reason is that the volume in the overdense regions 
is so small that their impact is negligible. In order to have slowdown, and after that 
acceleration, the fraction of the initial mass (or, equivalently, initial volume) which goes 
into the overdense structures should be larger. 

The peak model of structures. We have already mentioned that the distribution 
function (29) does not take into account the peak-in-a-peak problem, trough-in-a-trough 
problem or trough-in-a-peak problem. In particular, the last mentioned involves the 
extinction of voids by collapsing structures, which leads to the transfer of mass from 
underdense to overdense regions, increasing the volume fraction of the latter. In the 
excursion set formalism, the extinction of voids by overdense structures is crucial for 
getting a single void size to dominate at each epoch [73,127], as observed [139]. 

The peak model is based on an isolated view of structure formation, in which 
matter concentrations remain stationary, with mass accreting onto (or flowing away 
from) the fixed extremal points of the initial density field. Shear and tides are neglected. 
However, such effects can be important, and the sites of structure formation do not 
necessarily coincide with the initial density peaks (even underdense regions can collapse) 
[121,132,133,138]. 

The mass assignment of the peaks and troughs is likely to be important, and one 
should consider a realistic mass function instead of having all peaks and troughs contain 
the same amount of mass. Certainly, the treatment of merging and assignment of mass 
with the simple Gaussian smoothing is somewhat arbitrary, and the mass assignment 
should reflect the dynamics of structure formation. A realistic treatment would be 
expected to break the symmetry in mass between peaks and troughs. 

We have said that our model does not have new free parameters, and this is true 
in the sense that we have not adjusted anything to get a certain result. However, the 
smoothing threshold, which we set at o"o(t, R) = 1, could have been chosen to have some 
other value. The peak distribution (29) depends on 6/ao{t,R), and increasing ao{t,R) 
enhances the effect of structure. However, changing o"o(t, R) around unity does not make 
much difference. Setting ao{t,R) = 2 increases Ht by less than 10%, and the choice 
(Jo{t, R) = 0.5 brings Ht down by around 10%. The rate of change, i.e. the steepness of 
the Ht curve, increases with ao(t, R), and the deceleration parameter q changes by tens 
of percent when setting (Jo(t, R) to 2 or 0.5. The time when the expansion rate changes 
significantly remains about the same, and the qualitative behaviour is unchanged. 

Finally, the results depend strongly on the transfer function, and the BBKS 
approximation may not be sufficient near the important bend in the spectrum; this 
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is easy to correct in a numerical treatment. The results are not sensitive to the transfer 
function and power spectrum at small scales, where they are poorly known. The reason 
is that at early times when small structures form, their effect on the expansion rate is 
small. On the other hand, rapid changes in the power spectrum at large scales could 
have a large impact. Changing the primordial spectrum from scale-invariant to a power- 
law with a spectral index of n = 0.8 or n = 1.2 (with a pivot scale of 0.002 Mpc"^) 
changes Ht by only a few percent. More drastic changes in the power spectrum, such 
as a sudden transition, would be needed for a significant effect. 

A more realistic treatment of the effects discussed above will make it possible 
to say with more confidence whether overdense regions are prominent enough to give 
acceleration. If structures indeed lead to to accelerated expansion, the Buchert equations 
(7)-(9) dictate certain overall features that the universe must have, independent of the 
details of the structures. It is worth considering the compatibility of these general 
features with observations, as they might in principle rule out backreaction as an 
explanation for the observations. (See [29] for further comparison to observations.) 

5. Discussion 

5.1. Observational issues 

The age of the universe. Observationally, the age of the universe is not known with 
precision in a model-independent manner. There is an important constraint from the 
ages of globular clusters, which give the lower limit to > 11.2 Gyr at 95% C.L. and 
a best-fit age of to = 13.4 Gyr [140]. These results do not depend on the details of 
late-time cosmology, and should be valid as long as the expansion at redshifts 2; > 6 is 
close to the Einstein-de Sitter case. For the Hubble parameter, the most accurate 
model-independent determination is from SNe la observed with the Hubble Space 
Telescope [141]. The usually quoted value is h = 0.72 ± 0.08 [125] (in close agreement 
with h = 0.73 ± 0.06 found in [142]), while recent work with a different treatment of 
Cepheid metallicity gives h = 0.62 ± 0.05 [126] (all la limits). These estimates give a 
mean value of Hoto ~ 0.99 or i^o^o ~ 0.85, respectively, using to = 13.4 Gyr. Using the 
lower limit for to and the mean value for h gives Hoto > 0.83 or Hoto > 0.71, respectively. 
The treatment of Cepheids does not appear to be settled [141,143], but in any case, Hoto 
values in the lower range are not ruled out. In principle, Hoto can provide an important 
constraint on backreaction models. In particular, a definitive measurement of Ht > 1 
would imply that exotic matter is needed, because in a dust universe Ht < 1, regardless 
of the structures present (assuming that vorticity can be neglected) [129]. 

The values of Ht today in figure 4 are rather low for the BBKS transfer function, 
while the BD transfer function gives results more in agreement with the observations. 
The difference is mainly due to the fact that the change in the fraction of mass which is 
in structures is slower in the BBKS case. For the BBKS transfer function, only 30% of 
the mass is in structures at 15 Gyr, compared with 54% for the BD transfer function. 
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In both cases, initially (i.e. in the limit t 0) the fraction is 10%. 

The CMB peaks. As the exact expression (15) shows and figure 6 demonstrates, 
the average spatial curvature is generically large in a dust universe with significant 
dumpiness. This raises the question of compatibility of the backreaction explanation 
for the acceleration with CMB observations. However, the fact that the CMB data is in 
good agreement with FRW models that have no spatial curvature does not imply that 
the CMB rules out models with spatial curvature. Conclusions drawn about spatial 
curvature from the CMB are model- and prior- dependent. 

In physical terms, the CMB peak structure is mostly determined by the primordial 
spectrum of perturbations, by Q^h^ and Qi^h"^, and by the shift parameter, which is a 
measure of the distance to the last scattering surface [8,9]. The CMB anisotropies will 
be the mostly the same in models with identical values of these parameters, apart from 
the low multipoles (though see [144,145]). 

Even in the ACDM model with primordial perturbations given by a power law, the 
CMB alone does not constrain spatial curvature to be small [3], though combining the 
CMB with a measurement of Hubble parameter does provide a strong constraint. In a 
FRW model with a time-varying equation of state, at least ~ 0.2 is compatible with 
the correct CMB shift parameter (and fitting the SN la data and the A-parameter from 
baryon acoustic oscillations) [146]. 

In a backreaction model, there is no simple argument for obtaining the position of 
the CMB peaks. Both positively and negatively curved regions are present, and there 
may be more cancellation in the effect on the passage of light than would be expected 
based on a FRW model with the equivalent amount of spatial curvature. One has to redo 
the CMB analysis, considering the passage of light in a universe with realistic structures, 
as discussed in section 2.1. It might even be that the change in the passage of light 
together with the increase in Ht could explain the observations without acceleration, as 
discussed in [47]. 

In any case, the contribution of spatial curvature is regionally not negligible in the 
real universe. (The following discussion neglects effects due to breakdown of the dust 
approximation.) For example, consider an overdense dust structure. For a shell that is 
turning around from expansion to collapse, the left-hand side of (3) is zero. Neglecting 
vorticity (in the spherically symmetric case, it would be zero), the spatial curvature 
must be positive, and equal to the sum of the contributions of the energy density and 
shear, ^^^R = IGttGnP + 2cr^. For a stabilised structure, we have ^ = ^ = 0, so solving 
from (2) and (3) we obtain ^^''R = 127iG^p and u;^ = + cr^. Just as vorticity is 

needed to balance against the energy density and shear to have zero acceleration, zero 
expansion requires the spatial hypersurfaces to be positively curved. The energy density 
of stabilised structures can be much larger than the average energy density, so the same is 
true of their spatial curvature. (The growth in the relative spatial curvature is balanced 
by the fact that the volume occupied by stabilised regions shrinks by the same factor 
that their energy density increases, so the contribution to the average Hubble rate (8) 
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remains constant.) 

Inside voids, the energy density is mucli smaller than the background value, while 
the expansion rate is larger, so the contribution of the spatial curvature is again larger 
than that of energy density. For example, consider a void on an Einstein-de Sitter 
background. Neglecting shear and vorticity, we have I'-^-'-Rvoidl > ISttGnI^Kp), where 
6 is the real (not linear) density contrast of the void (typically 6 ~ —0.9 for observed 
voids [139]). 

In principle, observations of structures correlated with the CMB at different 
redshifts could lead to useful constraints on the evolution of the spatial curvature with 
redshift, as with the ISW effect [20], once the effect of non-linear structures on the 
passage of light is properly understood. 

Variance. A variance \fiQ\ > 0.2 is required to explain the observed acceleration [29]. 
The directional variation of the expansion rate was studied in [147] using SNe la 
(following the treatment of [125] rather than [126]). Typical variation was found at 
the 10-20% level, the maximum difference in the expansion rates being over 50%. Of 
course, angular variation is different from volume variance, and the systematics of SNe 
la are perhaps not completely understood, as discussed in section 1. For directional 
analysis of SNe la, see also [148,149]. 

In [150], the directional variation of the expansion rate inferred from SNe la in [125] 
was interpreted as evidence for backreaction. However, the estimate of the variance 
in [150] is unreliable (even apart from the applicability of linear perturbation theory). 
The magnitude of the the relevant term B = {di{di(pV'^(p)) — {di{dj(pdidj(p)) — |(V^v')^ 
was estimated by replacing spatial derivatives with the inverse of the averaging scale, 
and replacing the two powers of (f by the primordial amplitude of the power spectrum. 
However, spatial gradients are determined by the scale over which the perturbations 
vary significantly (as determined by the power spectrum and the transfer function), not 
by the averaging scale. For example, (5^) oc (V^v^V^v') is divergent for a scale-invariant 
spectrum (with no free-streaming cut-off), regardless of how large the averaging scale 
is, because there are perturbations on arbitrarily small scales. Conversely, with a free- 
streaming cut-off, the quantity is finite no matter how small the averaging scale is. 

Note that the variance in the expansion rate evaluated from Newtonian simulations 
[151] is large enough to provide acceleration, were it not balanced by shear, as noted 
in [29]. As we discuss in section 5.2 below, such a cancellation is not in present in general 
relativity. If the variance in simulations were negligible, it would be less plausible that 
it is large in the real relativistic case. However, the Newtonian variance is large, so in 
order to have acceleration, the shear simply has to be smaller than in the Newtonian 
case. 

It is sometimes argued that the dust shear would be of the order < 10~^°if^ 
from the isotropy of the CMB [152] (see also [153]). Any real non-linear dust structure 
violates this bound, so it is clear that the bound does not apply to the real universe. 
(The calculation also gives a limit of |Vp/p| < 10~^H for the spatial scale of density 
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variation, which is even more obviously violated.) However, if the shear really was 
negligible, then this would indicate that backreaction is important, as observations and 
simulations suggest a significant variance of the expansion rate [29]. 

Coldness of the Hubble flow. We have noted that there is no evidence for a 
matter component with negative pressure apart from the cosmological observations of 
accelerated expansion. In contrast to this, it has been claimed that the infiuence of 
vacuum energy can be seen in the local dynamics within the local 10 Mpc or so from 
the 'coldness' of the local Hubble fiow. 

The coldness refers to two distinct phenomena"*". First, the velocity dispersion in the 
local volume has been claimed to be anomalously low. Some authors find the dispersion 
within the nearest few Mpc to be ~ 40 km/s [154], while others find a value of 88 km/s 
lb 20 km/s [155], or over 100 km/s [156]. The second aspect is that the expansion rate 
measured within the nearest 10 Mpc is said to agree quite well with the global Hubble 
rate. Given that the matter distribution is locally quite clumpy, and does not become 
homogeneous until around 100 h~^Mpc, this seems somewhat surprising. It has been 
suggested that both observations are explained by vacuum energy. It can 'cool' the local 
expansion rate, and if vacuum energy dominates also the global dynamics, it is natural 
that the local and global expansion rates agree. 

The velocity dispersion argument was presented as evidence against a high density 
matter-dominated FRW model before the SN la observations supported accelerating 
expansion [157]. It has been argued that the observed low value of the velocity dispersion 
could be due to the known peculiarity of the local structures instead [158], or even that 
such a small dispersion is typical [159]. However, it seems that the typical velocity 
dispersion in simulations of the EdS model is indeed considerably higher, 300-700 
km/s [160], while ACDM simulations can reproduce the small velocity dispersion [161]. 
In simulations constrained to reproduce the large scale structure of the local universe, 
the velocity dispersion around regions similar to the Local Group of galaxies has been 
found to be as cold in the open CDM model as in the ACDM model, given the same 
matter density [162]. This still leaves the question of how common such regions are. 
However, the lower range of the values 150-300 km/s found in unconstrained simulations 
of the open CDM model [160] is not very dissonant with an observational value of over 
100 km/s. It does not seem unreasonable that the effect of spatial curvature would be 
even stronger in a clumpy model, where the density parameter of spatial curvature can 
be larger. 

Unlike for the velocity dispersion, the probability of being located in a region where 
the local expansion rate agrees with the global value has not been studied quantitatively. 
At any rate, as long as the measurement of the global and local Hubble parameters is 
uncertain, it seems difficult to draw strong conclusions. The notably different values 
h = 0.62 and h = 0.72 have both been used as evidence for the vacuum energy origin of 

+ The linearity of the nearby mean flow with distance is sometimes mentioned as a third aspect. 
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the local expansion rate [126, 163, 164]. 

It has also been argued that the influence of vacuum energy has been directly 
measured in the motions of nearby galaxies [163, 164]. The idea is that for a spherically 
symmetric system with a large central mass, the repulsive gravity of vacuum energy 
dominates beyond a certain radius, estimated as 1-2 Mpc for the Local Group. However, 
the assumption of a spherically symmetric field generated by a point mass used in 
the analysis does not hold very well, and the observations are also consistent with 
domination by negative curvature [162]. 

5.2. Non- Newtonian effects 

Spatial curvature and Newtonian gravity. It is not clear how much the overdense regions 
slow down the expansion rate, but from the model we have discussed, it seems difficult 
to avoid the conclusion that underdense regions cause Ht to increase significantly. The 
physical interpretation seems straightforward: the fraction of space in faster expanding 
regions grows. However, there is a subtlety involved. The evolution we took for the 
individual regions is purely Newtonian, and the peak statistics do not involve general 
relativity. Had we formulated the problem in Newtonian gravity instead of general 
relativity, the answer would have been the same. However, in Newtonian gravity, the 
backreaction variable Q is a boundary term [76], so backreaction vanishes for periodic 
boundary conditions. (In particular, this is true for Newtonian simulations.) 

While the universe presumably does not have periodic boundary conditions at the 
visual horizon, the result implies that backreaction in Newtonian gravity also vanishes 
for a statistically homogeneous and isotropic system. This can be seen in two ways. 
First, consider a volume which has periodic boundary conditions on a scale much larger 
than the homogeneity scale. By statistical homogeneity and isotropy, the value of Q 
evaluated within each homogeneity scale sized box in the volume is the same as the 
overall value, which is zero due to the boundary conditions. Because the local physics 
is independent of the boundary conditions on very large scales, this result should hold 
even if the boundary conditions are not periodic*. A perhaps more physical argument, 
explained in [128], is that a boundary term can be viewed as a fiux going from one 
region to another, and in a statistically homogeneous and isotropic region, there should 
be an equal fiux in and out. 

In Newtonian gravity, the average evolution of a statistically homogeneous and 
isotropic dust space does always follow the FRW equations (i.e. the Buchert equations 
with Q = 0). However, this is not true in general relativity. The Newtonian theory 
constraint that variance and shear in Q in (10) cancel up to a boundary term is not 
present in general relativity [77]. Viewed equivalently in terms of the spatial curvature 
via the integrability condition (11), this difference is a refiection of the absence of spatial 
curvature in Newtonian gravity. 

In Newtonian space and time, the geometry of the spatial hypersurfaces is fixed, 

* I am grateful to Christof Wetterich for this argument. 
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so the spatial curvature tensor and its trace, the spatial curvature scalar ^^^R, do not 
exist [99]. Therefore the Newtonian equations of motion do not include an equation 
for the expansion rate such as the Hamiltonian constraint (3) independent of the 
Raychaudhuri equation. The analogue of the relativistic Hubble equation (3) emerges 
only as the first integral of the Raychaudhuri equation. 

In Newtonian gravity, the backreaction variable Q is given by a boundary term 
which is in general non-zero. However, for an isolated system or one that is statistically 
homogeneous and isotropic, we have Q = 0. Then the the first integral of the Newtonian 
equivalent of the average Raychaudhuri equation (7) gives an equation like the average 
Hamiltonian constraint (8), but with a term proportional to in place of {^^^R), as 
(14) shows. The term multiplied by can be interpreted as the conserved energy 
of the Newtonian system. 

There is no such conserved quantity in the relativistic case. The average spatial 
curvature can evolve non-trivially, unlike the total energy of an isolated Newtonian 
system. The average spatial curvature is constrained by the integrability condition 
(11) for the average Hamiltonian constraint and the average Raychaudhuri equation, 
instead of being proportional to a~^. (In the FRW case, the relativistic spatial curvature 
term and the Newtonian energy term are mathematically identical, only the physical 
interpretation is different.) 

In our model, the evolution of overdense and underdense regions is uncorrelated. We 
do not have a Newtonian constraint on their behaviour, and the overall average spatial 
curvature can evolve non-trivially, as in the general relativistic case. A Newtonian 
calculation for a statistically homogeneous and isotropic system would have to include 
the global constraint that the underdense and overdense regions add up in such a manner 
that the energy term (corresponding to the spatial curvature) is proportional to a~^. 
(In a consistent treatment with a continuous distribution of matter, instead of isolated 
regions, this would follow automatically.) 

There is a similar situation in the case of spherical symmetry. As noted in 
section 4.1, the backreaction variable Q vanishes for spherical symmetry in Newtonian 
gravity [26,78], resulting in the spherical collapse/expansion model. In general relativity, 
spherical symmetry does not imply perfect cancellation between variance and shear, 
and the average expansion of a spherically symmetric dust solution can even accelerate 
[30,113,114]. 

The Newtonian limit of general relativity. In order for backreaction to be able to 
account for the observed acceleration, non-Newtonian aspects of gravity should be 
important already at the homogeneity scale of around 100 h~^Mpc. It is sometimes 
argued that general relativistic effects in the present-day universe can only be important 
on super- Hubble scales or near neutron stars and black holes. However, this conclusion is 
based on linearly perturbed FRW or Minkowski universes or the Schwarzschild solution, 
and the general situation is not that simple. (The subtleties of the Newtonian limit of 
near-FRW cosmological models have been discussed in [165,166]; see also [167] regarding 
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the magnetic part of the Weyl tensor.) 

The Einstein equation has ten components, with four constraints, whereas in 
Newtonian gravity there is only one equation, the Poisson equation. There are two 
aspects to the additional equations. Some of them correspond to degrees of freedom 
which do not exist in Newtonian gravity (such as gravity waves and spatial curvature), 
while others provide new constraints on the existing degrees of freedom. The low- 
velocity, weak-field limit of general relativity is not Newtonian gravity, but Newtonian 
gravity with extra degrees of freedom and additional constraints. 

General relativity reduces to Newtonian gravity in the limit of infinite speed of 
light. However, this limit is not smooth. Taking the speed of light to infinity removes 
some equations completely and makes Newtonian gravity qualitatively different from 
general relativity, where the speed of light is finite. Relativistic effects do not require 
velocities near the speed of light, only a finite speed of hght. 

For example, there exist solutions of Newtonian gravity which are not the limit 
of any general relativity solution, due to additional constraints in the relativistic case. 
There are Newtonian expanding dust solutions which have zero shear but non-zero 
vorticity. However, in general relativity, non-zero vorticity implies non-zero shear for 
expanding dust [99,106]. This does not require large velocities or strong gravitational 
fields; see [168] for a clear comparison of the general relativistic and Newtonian cases. 

In the case of backreaction, it is the presence of new degrees of freedom, namely 
spatial curvature, rather than extra constraints which is important. In perturbation 
theory around FRW universes, the leading terms in the backreaction variable Q are 
Newtonian, so they are total derivatives and do not contribute in a statistically 
homogeneous and isotropic system [25,81,128]. The non- Newtonian terms do not appear 
earlier than fourth order in perturbation theory. Their numerical coefficients have not 
been evaluated (it may be possible to do so using third order perturbation theory, instead 
of having to calculate to fourth order [169]), but their form is known, and velocities near 
the speed of light are not required for them to be important. 

Simulations and backreaction. If non-Newtonian aspects are important already at the 
homogeneity scale, one would expect them to show up also in quantities other than 
the expansion rate. We can ask whether there is any room for such effects, given the 
comparison of simulations of structure formation (which are completely NeAvtonian) 
with observations. In fact, while simulations reproduce many features of observations, 
there are some notable differences. 

In [170] it was found that the homogeneity scale evaluated from simulations was 
only of the order 10 h~^Mpc, an order of magnitude smaller than the observed 70-100 
h~^Mpc [89,90]. (The box size of the simulations studied in [170] was 141 h"^Mpc, 
so the correct homogeneity scale could not have been reliably recovered. However, the 
homogeneity scale that was found is much smaller than the box size, so it is not likely 
to be a finite size artifact.) 

According to [171], the number of largest structures in the Millennium Simulation 
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is underproduced compared to observations, by a factor of 10 for the most luminous 
super clusters. As for large underdense structures, the voids produced in simulations are 
not as empty as those observed. This 'void phenomenon' has even been called a crisis 
of the ACDM model [172]. (It has recently been argued that simulations coupled with 
semianalytical galaxy formation in fact agree with the void observations [173].) 

These discrepancies may have resolutions which have nothing to do with non- 
Newtonian physics or backreaction. The statistics for the largest structures are probing 
the tail of the distribution, which may be affected by non-linear corrections to the 
evolution of the power spectrum [97] or transients from initial conditions [174]. The 
magnitude of such effects has been found at the 10-30% level, and the discrepancy is an 
order of magnitude, but one should not be confident that there are no unaccounted for 
effects in simulations, for example related to discretisation [95]. The apparent emptiness 
of voids, in turn, might be due to bias or subtleties in galaxy formation and the definition 
of voids, rather than spatial curvature [175,176]. 

Nevertheless, the point is that we do not have confirmation that Newtonian gravity 
works well near the homogeneity scale, since the Newtonian results differ in some 
important respects from the observations. 

5.3. Conclusion 

Summary. We have studied the effect of structure formation on the expansion rate 
in a statistically homogeneous and isotropic space. We first reviewed studies of the 
propagation of light in a space with non-linear structures, and discussed some qualitative 
issues related to the average expansion rate. We then calculated the average expansion 
rate in a model where the number density of structures is given by the peak model of 
structure formation with cold dark matter, and the individual structures are described 
with the spherical collapse model and its underdense equivalent. In the calculation, there 
are no adjustable free parameters in addition to those related to structure formation in 
a spatially flat matter-dominated FRW universe. 

We find that the expansion rate increases relative to the FRW value at a time of 
tens of billions years, about the observed acceleration era, possibly offering a solution 
to the coincidence problem. The timescale has its origin in the change of slope of the 
CDM transfer function around the matter-radiation equality scale fceq. This leads to a 
preferred time around lO^teq ~ 10^° years, when the volume occupied by structures and 
their size relative to the visual horizon saturate, and the impact of structures on the 
expansion rate becomes large. 

However, while Ht increases in the model, it does not rise sufficiently rapidly to 
correspond to acceleration. The expansion rate increases because the relative volume 
of the faster expanding regions rises, as quantified by the Buchert equations. In order 
to have acceleration, it is likely that the average expansion rate would first have to be 
damped by slower expanding overdense regions before it is increased by shedding their 
contribution when voids become dominant. In the present model, this does not happen. 
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because the volume occupied by the overdense regions is rather small. This is partly due 
to the fact that we have equal amounts of mass in overdense and underdense regions. 
Taking into account mass flow to the overdense regions would increase their impact. 

Improving the treatment of the statistics of the peaks and the model used for 
individual structures would lead to a more realistic estimate of the expansion rate. The 
propagation of light in a universe with realistic, evolving structures also has to be further 
studied to see how the average expansion rate is related to observations. 

If a realistic backreaction model does turn out to provide acceleration in agreement 
with the observations, one could say that it unifies three historically popular FRW 
models. There is only matter present and the initial value of the spatial curvature at 
the background level is zero as in the standard CDM model, the universe has negative 
curvature as in the open CDM model, and the expansion accelerates as in the ACDM 
model. Acceleration due to structure formation would probably be a transient phase on 
the way to negative curvature voids dominating the expansion of the universe, though the 
behaviour at very late times depends on how the universe is structured on scales which 
are currently far beyond the horizon, of which we have no theoretical or observational 
understanding. 
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